Tutorial Functions average distance and average density

On this page

• 1 Tutorial to calculate average distance and point density statistics
• 2 Overview of Steps
• 3 Loading packages
• 4 Setting up working directory and path to wrapper functions
• 5 Calculating average distance to rivers for each polygon ID
• 6 Loading inputs
• 7 Calculating average distance to rivers for each ID polygon in target
• 8 Calculating average density to point-buildings for each ID polygon in target
• 9 Using Buildings dataset as shape_input
• 10 Using points rather than non-intersecting buffers

1 Tutorial to calculate average distance and point density statistics

This guide provides a step-by-step approach to compute the average distance and point density statistics for a target set of polygons/points coordinates in R. It uses own functions from climatic4economist package: ave_dist_2shp() and ave_point_density_2shp(). The inputs are shapefile files loaded as sf_dataframe R objects and the returning outputs are dataframe R objects.

2 Overview of Steps

We will go through the following steps:

1. Load the data
2. Prepare the data
3. Compute the average distance to rivers for each ID (village, adm-div, etc.)
4. Compute the average point density of buildings for each ID (village, adm-div, etc.)

3 Loading packages

library("dplyr")
library("raster")
library("sf")
library("data.table")
library("datawizard")
library("zoo")
library("doParallel")
library("foreach")
library("parallel")
library("ncdf4")
library("SPEI")
library("ggplot2")
library("foreign")
library("tidyr")
library("terra")
library("stars")
library("yarrr")
library("wesanderson")
#library("ClimateOperators")
library("rlist")
library("writexl")
library("readxl")
library("stringr")
library("tidyr")
library("miceadds")
library("stargazer") #Display table on latex format
library("foreign")
library("dplyr") #data management
library("tidyverse") #data management
library("tibble") # tibbles
library("caret")
library("naniar")
library("mapproj")
library("ggmap")
library("RgoogleMaps")
library("mapview")
library("leaflet")
library("haven")
library("rnaturalearth")
library("reshape2")
library("sjlabelled")
library("RColorBrewer")
library("climatic4economist")
devtools::document() 

All packages loaded

4 Setting up working directory and path to wrapper functions

Working directory set up as default with folder’s location.

# Otherwise use
#setwd("C:/Users/X/OneDrive - Food and Agriculture Organization/PROJECT11_GEOINDICATORS_REPOSITORY")

rm(list = ls()) # remove existing objects

# load wrapper functions
# source(file.path(#"web_tutorial",
#                  "functions.R"))

path_to_project <- file.path("..", "..",
                        "data", "data_project", "benin",
                        "PROJECT11_GEOINDICATORS_REPOSITORY")

5 Calculating average distance to rivers for each polygon ID

6 Loading inputs

All input have to be set up as EPSG4326 or WGS84 as coordinate system

6.0.1 Shapefile extension for the target country (Benin)

This uses GADM database (see GADM)

path_gadm3 <- file.path(path_to_project,
                        "INPUT",
                        "Benin_shp","Benin-ADMlevels","Shapefile",
                        "BEN_shp","gadm41_BEN_3.shp")

gadm3 <- st_as_sf(st_read(path_gadm3)) |>
  dplyr::select(-NL_NAME_1, -NL_NAME_2, -NL_NAME_3, 
                -VARNAME_3, -CC_3, -HASC_3,-TYPE_3,
                -ENGTYPE_3,-COUNTRY)

gadm3 <- st_transform(gadm3, 4326) #defining as EPSG4326 or WGS84
gadm3

Simple feature collection with 546 features and 7 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 0.774345 ymin: 6.23491 xmax: 3.851701 ymax: 12.41835
Geodetic CRS:  WGS 84
First 10 features:
          GID_3 GID_0   GID_1  NAME_1     GID_2    NAME_2     NAME_3
1   BEN.1.1.1_1   BEN BEN.1_1 Alibori BEN.1.1_1 Banikoara  Banikoara
2   BEN.1.1.2_1   BEN BEN.1_1 Alibori BEN.1.1_1 Banikoara   Founougo
3   BEN.1.1.3_1   BEN BEN.1_1 Alibori BEN.1.1_1 Banikoara   Gomparou
4   BEN.1.1.4_1   BEN BEN.1_1 Alibori BEN.1.1_1 Banikoara    Goumori
5   BEN.1.1.5_1   BEN BEN.1_1 Alibori BEN.1.1_1 Banikoara      Kokey
6   BEN.1.1.6_1   BEN BEN.1_1 Alibori BEN.1.1_1 Banikoara Kokiborou*
7   BEN.1.1.7_1   BEN BEN.1_1 Alibori BEN.1.1_1 Banikoara      Ounet
8   BEN.1.1.8_1   BEN BEN.1_1 Alibori BEN.1.1_1 Banikoara Somp?r?kou
9   BEN.1.1.9_1   BEN BEN.1_1 Alibori BEN.1.1_1 Banikoara     Soroko
10 BEN.1.1.10_1   BEN BEN.1_1 Alibori BEN.1.1_1 Banikoara      Toura
                         geometry
1  MULTIPOLYGON (((2.415323 11...
2  MULTIPOLYGON (((2.346528 11...
3  MULTIPOLYGON (((2.464541 11...
4  MULTIPOLYGON (((2.070032 11...
5  MULTIPOLYGON (((2.802918 11...
6  MULTIPOLYGON (((2.26218 11....
7  MULTIPOLYGON (((2.405739 11...
8  MULTIPOLYGON (((2.463629 11...
9  MULTIPOLYGON (((2.338276 11...
10 MULTIPOLYGON (((2.317923 11...

6.0.2 Shapefile target country (Benin)

This file is the shapefile with the polygon’s structure. Each ID is enumerated in column id_gns. The original file with X, Y villages’ coordinates can be found here Villages Benin. The file was modified and non-intersecting buffers (using Voronoi’s structure) of 2km were constructed using the villages’ coordinates.

path_polygons_gns <- file.path(path_to_project,
                        "INPUT","QGIS_map_voronoi",
                        "villages_voronoibuffer_adm4_gns.shp")

polygons_gns <- st_as_sf(st_read(path_polygons_gns))

polygons_gns <- st_transform(polygons_gns, 4326) #defining as EPSG4326 or WGS84
polygons_gns

Simple feature collection with 3599 features and 18 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 0.7645527 ymin: 6.222888 xmax: 3.839671 ymax: 12.42317
Geodetic CRS:  WGS 84
First 10 features:
   ID gid_0   gid_1     gid_2       gid_3  name_1    name_2    name_3
1   1   BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
2   1   BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
3   1   BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
4   1   BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
5   1   BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
6   1   BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
7   1   BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
8   1   BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
9   1   BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
10  1   BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
           type_3 engtype_3 id_gns  full_name full_nm_nd sort_name desig_cd
1  Arrondissement   Borough   2058       Gano       Gano      GANO      PPL
2  Arrondissement   Borough   1426     Kokabo     Kokabo    KOKABO      PPL
3  Arrondissement   Borough    973      Wagui      Wagui     WAGUI      PPL
4  Arrondissement   Borough   3392   Atabénou   Atabenou  ATABENOU      PPL
5  Arrondissement   Borough   2726      Komon      Komon     KOMON      PPL
6  Arrondissement   Borough   2263      Dérou      Derou     DEROU      PPL
7  Arrondissement   Borough   3386  Banikoara  Banikoara BANIKOARA      PPL
8  Arrondissement   Borough    453 Toké-Banta Toke-Banta TOKEBANTA      PPL
9  Arrondissement   Borough   1873   Gomparou   Gomparou  GOMPAROU      PPL
10 Arrondissement   Borough    581    Sonworé    Sonwore   SONWORE      PPL
     adm1    lat_y    lon_x                       geometry
1  BJ-000 11.24685 2.414291 MULTIPOLYGON (((2.423529 11...
2  BJ-000 11.28007 2.427724 MULTIPOLYGON (((2.434365 11...
3  BJ-000 11.28146 2.400951 MULTIPOLYGON (((2.392831 11...
4  BJ-000 11.28761 2.381083 MULTIPOLYGON (((2.392831 11...
5  BJ-000 11.28951 2.423543 MULTIPOLYGON (((2.426221 11...
6  BJ-000 11.29832 2.401679 MULTIPOLYGON (((2.413442 11...
7   BJ-AL 11.29845 2.438561 MULTIPOLYGON (((2.447332 11...
8  BJ-000 11.30770 2.414829 MULTIPOLYGON (((2.430931 11...
9   BJ-AL 11.33134 2.441141 MULTIPOLYGON (((2.437974 11...
10 BJ-000 11.35827 2.422149 MULTIPOLYGON (((2.440276 11...

6.0.3 Shapefile shape_input

We will use the shapefile of African Rivers from FAO as our shape_input (see FAO rivers1 and here FAO rivers2)

path_rivers_fao <- file.path(path_to_project, "INPUT","FAO_GIS","rivers_africa","rivers_africa_37333.shp")

rivers_fao <- st_as_sf(st_read(path_rivers_fao))

rivers_fao <- st_transform(rivers_fao, 4326) #defining as EPSG4326 or WGS84
rivers_fao

Simple feature collection with 185730 features and 18 fields
Geometry type: MULTILINESTRING
Dimension:     XY
Bounding box:  xmin: -17.23958 ymin: -34.76458 xmax: 54.35 ymax: 37.21875
Geodetic CRS:  WGS 84
First 10 features:
   FID_af_str ARCID FROM_NODE TO_NODE FID_sub_ba SUB_BAS MAJ_BAS
1           0     1         2       1      32890  201591    7020
2           1     2         5       4      32890  201591    7020
3           2     3         6       7      32890  201591    7020
4           3     4         7       3      32890  201591    7020
5           4     5         8       6      32890  201591    7020
6           5     6        10       8      32890  201591    7020
7           6     7        11       9      32890  201591    7020
8           7     8        13       8      32890  201591    7020
9           8     9        14       6      32890  201591    7020
10          9    10        15       7      32890  201591    7020
                    MAJ_NAME            SUB_NAME MAJ_AREA LEGEND SUBBAS_ID
1  Mediterranean South Coast Algerian east coast   558292     20   7201591
2  Mediterranean South Coast Algerian east coast   558292     20   7201591
3  Mediterranean South Coast Algerian east coast   558292     20   7201591
4  Mediterranean South Coast Algerian east coast   558292     20   7201591
5  Mediterranean South Coast Algerian east coast   558292     20   7201591
6  Mediterranean South Coast Algerian east coast   558292     20   7201591
7  Mediterranean South Coast Algerian east coast   558292     20   7201591
8  Mediterranean South Coast Algerian east coast   558292     20   7201591
9  Mediterranean South Coast Algerian east coast   558292     20   7201591
10 Mediterranean South Coast Algerian east coast   558292     20   7201591
   TOBAS_ID Strahler A_Strahler RASTERVA_2 RASTERVA_1 Regime
1      -999        1          8        609       5202      P
2      -999        1          8        460       2537      P
3      -999        2          7        518      97325      P
4      -999        3          6        478     218144      P
5      -999        2          7        518      41565      P
6      -999        1          8        518       8568      P
7      -999        1          8        413       9945      P
8      -999        1          8        638       5032      P
9      -999        1          8        648       6162      P
10     -999        2          7        478     106369      P
                         geometry
1  MULTILINESTRING ((9.220833 ...
2  MULTILINESTRING ((9.945833 ...
3  MULTILINESTRING ((9.65625 3...
4  MULTILINESTRING ((9.735417 ...
5  MULTILINESTRING ((9.647917 ...
6  MULTILINESTRING ((9.614583 ...
7  MULTILINESTRING ((10.075 37...
8  MULTILINESTRING ((9.479167 ...
9  MULTILINESTRING ((9.2375 37...
10 MULTILINESTRING ((9.68125 3...

As the shapefile of rivers corresponds to all Africa, we will crop it to country’s area using the extension from the object gamd, which has the extension of Benin. It will be enlarged a bit in case to guarantee it does contain all country’s borders. Notice also that for using the function raster::crop, we need to define the object as an Spatial one, do the cropping and define it again as sf_dataframe, for which, we need to define it again with the same coordinate system WGS84 (otherwise our function ave_point_density_2shp will trew an error).

ext_benin<-raster::extent(gadm3)
ext_benin@xmin<-ext_benin@xmin-0.5
ext_benin@xmax<-ext_benin@xmax+0.5
ext_benin@ymin<-ext_benin@ymin-0.5
ext_benin@ymax<-ext_benin@ymax+0.5

#cropping
rivers_fao <- st_as_sf(raster::crop(as(rivers_fao, "Spatial"),ext_benin))

rivers_fao <- st_transform(rivers_fao, 4326) #defining as EPSG4326 or WGS84 again as I only defined as sf dataframe before
rivers_fao

Simple feature collection with 1796 features and 18 fields
Geometry type: MULTILINESTRING
Dimension:     XY
Bounding box:  xmin: 0.274345 ymin: 5.78125 xmax: 4.351701 ymax: 12.91835
Geodetic CRS:  WGS 84
First 10 features:
   FID_af_str ARCID FROM_NODE TO_NODE FID_sub_ba SUB_BAS MAJ_BAS MAJ_NAME
1       85827 85828     89460   89298      32706   20872    7002    Niger
2       85828 85829     89461   89407      32705   20871    7002    Niger
3       85829 85830     89348   89462      32705   20871    7002    Niger
4       85871 85872     89334   89502      32706   20872    7002    Niger
5       85872 85873     89502   89460      32706   20872    7002    Niger
6       85918 85919     89264   89462      32705   20871    7002    Niger
7       85930 85931     89569   89381      32706   20872    7002    Niger
8       85931 85932     89070   89571      32692   20843    7002    Niger
9       85979 85980     88722   89617      32701   20854    7002    Niger
10      85980 85981     89553   89617      32701   20854    7002    Niger
        SUB_NAME MAJ_AREA LEGEND SUBBAS_ID TOBAS_ID Strahler A_Strahler
1           Faga  2136941      2   7020872  7020871        4          5
2        Niger 9  2136941      2   7020871  7020859        1          8
3        Niger 9  2136941      2   7020871  7020859        6          3
4           Faga  2136941      2   7020872  7020871        1          8
5           Faga  2136941      2   7020872  7020871        4          5
6        Niger 9  2136941      2   7020871  7020859        4          5
7           Faga  2136941      2   7020872  7020871        1          8
8       Sokoto 1  2136941      2   7020843  7020841        1          8
9  Dallol Maouri  2136941      2   7020854  7020853        2          7
10 Dallol Maouri  2136941      2   7020854  7020853        1          8
   RASTERVA_2 RASTERVA_1 Regime                       geometry
1         182     907940      P MULTILINESTRING ((0.5229167...
2         421      10060      I MULTILINESTRING ((2.102083 ...
3          55   66217700      P MULTILINESTRING ((2.376455 ...
4         667       5300      I MULTILINESTRING ((0.4721039...
5         674     894820      P MULTILINESTRING ((0.4979167...
6         158    1462580      P MULTILINESTRING ((2.256649 ...
7         262       5870      I MULTILINESTRING ((0.76875 1...
8         280       6090      I MULTILINESTRING ((4.34375 1...
9         705     178340      I MULTILINESTRING ((3.49375 1...
10        586       5100      I MULTILINESTRING ((3.514583 ...

7 Calculating average distance to rivers for each ID polygon in target

This wrapper function calculates the average distance (in m) of each ID element in target to the object shape_input. It uses distanceFromPoints which works better, faster than st_distance function. For the moment, the function works better for UTM coordinate system.

country_ext defines the whole extension of the country, including country borders. Here, we use gadm3 is our country’s extension for Benin.

target is the shapefile (adm-div such as communes, villages or HH coordinates) for which the average distance will be calculated. We use polygons_gns as our target polygon for which we will calculate the distance.

shape_input is the shapefile for which a raster of distances will be calculated (eg. roads, rivers, etc.). In our case, we use rivers_fao as the object from which the distances will be calculated.

UTM corresponds to the country’s location in the UTM coordinate system. As for calculating distances it is more convenient to use meters instead of degrees, and the coordinate system needs to be transformed into the UTM of the country, which is more precise around the equator. We define UTM=31, which corresponds to Benin. Otherwise, it uses default at the global level (3857) (under development).

res_raster is the resolution in meters for the raster of distances that will be created. A higher resolution increase the time of processing. We define a resolution of 10000 which corresponds to 10 km. Internally, it will create a raster of distances of 10 km, and each pixel size of 10kmX10km.

name_dist name of the average distance variable constructed. It will be created as ave_dist_NAME with name_dist=“NAME” giving the name of the distance variable that will be created. In our case, ave_dist_rivers.

ave_dist_rivers <- ave_dist_2shp(country_ext=gadm3,target=polygons_gns,shape_input=rivers_fao,UTM=31,res_raster=10000,name_dist="rivers")

[1] "UTM is: 31"

ave_dist_rivers[1:10, ] #print first 10 elements

   ID gid_0   gid_1     gid_2       gid_3  name_1    name_2    name_3
1   1   BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
2   1   BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
3   1   BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
4   1   BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
5   1   BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
6   1   BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
7   1   BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
8   1   BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
9   1   BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
10  1   BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
           type_3 engtype_3 id_gns  full_name full_nm_nd sort_name desig_cd
1  Arrondissement   Borough   2058       Gano       Gano      GANO      PPL
2  Arrondissement   Borough   1426     Kokabo     Kokabo    KOKABO      PPL
3  Arrondissement   Borough    973      Wagui      Wagui     WAGUI      PPL
4  Arrondissement   Borough   3392   Atabénou   Atabenou  ATABENOU      PPL
5  Arrondissement   Borough   2726      Komon      Komon     KOMON      PPL
6  Arrondissement   Borough   2263      Dérou      Derou     DEROU      PPL
7  Arrondissement   Borough   3386  Banikoara  Banikoara BANIKOARA      PPL
8  Arrondissement   Borough    453 Toké-Banta Toke-Banta TOKEBANTA      PPL
9  Arrondissement   Borough   1873   Gomparou   Gomparou  GOMPAROU      PPL
10 Arrondissement   Borough    581    Sonworé    Sonwore   SONWORE      PPL
     adm1    lat_y    lon_x ave_dist_rivers
1  BJ-000 11.24685 2.414291        4885.229
2  BJ-000 11.28007 2.427724        4885.229
3  BJ-000 11.28146 2.400951        4885.229
4  BJ-000 11.28761 2.381083        4885.229
5  BJ-000 11.28951 2.423543        4885.229
6  BJ-000 11.29832 2.401679        4885.229
7   BJ-AL 11.29845 2.438561        5402.293
8  BJ-000 11.30770 2.414829        4885.229
9   BJ-AL 11.33134 2.441141        6628.781
10 BJ-000 11.35827 2.422149        6761.864

8 Calculating average density to point-buildings for each ID polygon in target

This wrapper function calculates the average density of points for each ID element of targetto the object shape_input .

9 Using Buildings dataset as shape_input

While some of the objects are already defined as in the previous function (gadm3 and polygons_gns), we will use as shape_input a shapefile that contains the points with buildings located in the country (see OCHA Datahum). The file comes in zip format as it is opened, and transformed as an sf_dataframe

path2buildings <- file.path(path_to_project,
                        "INPUT", "Benin_OCRA","Buildings","hotosm_ben_buildings_polygons_shp.zip") #using zip

# Extract the ZIP file and Find the file among the extracted files. There is only one
output_directory_buildings <- tempdir()  # Using a temporary directory
unzip(path2buildings, exdir = output_directory_buildings)
extracted_files <- list.files(output_directory_buildings, full.names = TRUE)
buildings_file_path <- extracted_files[grep("\\.shp$", extracted_files)]

#reading file
buildings_ocra <- st_as_sf(st_read(buildings_file_path))

buildings_ocra <- st_transform(buildings_ocra, 4326) #defining as EPSG4326 or WGS84
buildings_ocra

Simple feature collection with 813093 features and 11 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 0.7848944 ymin: 6.236586 xmax: 3.836711 ymax: 12.40539
Geodetic CRS:  WGS 84
First 10 features:
      osm_id addrcity office                      name buildingle addrhousen
1   81452793     <NA>   <NA>                Tri Postal       <NA>       <NA>
2   81766299     <NA>   <NA>                      <NA>          8       <NA>
3  116969223     <NA>   <NA>                Site Pompe       <NA>       <NA>
4  116969224     <NA>   <NA>             Station SONAB       <NA>       <NA>
5  116969225     <NA>   <NA>                Site Pompe       <NA>       <NA>
6  116985473     <NA>   <NA>                      SBEE       <NA>       <NA>
7  116985474     <NA>   <NA> Hôtel de Ville de Bassila       <NA>       <NA>
8  116985476     <NA>   <NA> Hôtel de Ville de Bassila       <NA>       <NA>
9  167008127     <NA>   <NA>                      <NA>       <NA>       <NA>
10 167008128     <NA>   <NA>                      <NA>       <NA>       <NA>
   addrfull building addrstreet source buildingma
1      <NA>      yes       <NA>   <NA>       <NA>
2      <NA>      yes       <NA>   <NA>       <NA>
3      <NA>      yes       <NA>   <NA>       <NA>
4      <NA>      yes       <NA>   <NA>       <NA>
5      <NA>      yes       <NA>   <NA>       <NA>
6      <NA>      yes       <NA>   <NA>       <NA>
7      <NA>      yes       <NA>   <NA>       <NA>
8      <NA>      yes       <NA>   <NA>       <NA>
9      <NA>      yes       <NA>   <NA>       <NA>
10     <NA>      yes       <NA>   <NA>       <NA>
                         geometry
1  MULTIPOLYGON (((2.386848 6....
2  MULTIPOLYGON (((2.386307 6....
3  MULTIPOLYGON (((1.666625 9....
4  MULTIPOLYGON (((1.667526 9....
5  MULTIPOLYGON (((1.6682 9.00...
6  MULTIPOLYGON (((1.667723 8....
7  MULTIPOLYGON (((1.667386 8....
8  MULTIPOLYGON (((1.667395 8....
9  MULTIPOLYGON (((2.340467 6....
10 MULTIPOLYGON (((2.340359 6....

Below we describe the inputs for the function:

country_ext defines the whole extension of the country, including country borders. Here, we use gadm3 is our country’s extension for Benin.

target is the shapefile (adm-div such as communes, villages or HH coordinates) for which the average point density will be calculated. We use polygons_gns as our target polygon for which we will calculate the point density.

shape_input is the shapefile from which to calculate points density (eg. buildings, roads, rivers, etc.) In our case, it will be our building’s object buildings_ocra.

res_raster is the resolution in degrees for the raster of distances that will be created and applied to `shape_input’ (0.1 could be equivalent to 10km or 0.1*100000=10000m around the Equator). A higher resolution increases the time of processing.

name_dist name of the average density variable constructed. It will be created as ave_den_NAME with name_density=“NAME In our case, ave_den_allbuildings.

ave_den_allbuildings <-ave_point_density_2shp(country_ext=gadm3,
                                              target=polygons_gns,
                                              shape_input=buildings_ocra,
                                              res_raster=0.1,
                                              name_density="allbuildings")
ave_den_allbuildings[1:10, ] #print first 10 elements

   ID gid_0   gid_1     gid_2       gid_3  name_1    name_2    name_3
1   1   BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
2   1   BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
3   1   BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
4   1   BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
5   1   BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
6   1   BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
7   1   BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
8   1   BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
9   1   BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
10  1   BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
           type_3 engtype_3 id_gns  full_name full_nm_nd sort_name desig_cd
1  Arrondissement   Borough   2058       Gano       Gano      GANO      PPL
2  Arrondissement   Borough   1426     Kokabo     Kokabo    KOKABO      PPL
3  Arrondissement   Borough    973      Wagui      Wagui     WAGUI      PPL
4  Arrondissement   Borough   3392   Atabénou   Atabenou  ATABENOU      PPL
5  Arrondissement   Borough   2726      Komon      Komon     KOMON      PPL
6  Arrondissement   Borough   2263      Dérou      Derou     DEROU      PPL
7  Arrondissement   Borough   3386  Banikoara  Banikoara BANIKOARA      PPL
8  Arrondissement   Borough    453 Toké-Banta Toke-Banta TOKEBANTA      PPL
9  Arrondissement   Borough   1873   Gomparou   Gomparou  GOMPAROU      PPL
10 Arrondissement   Borough    581    Sonworé    Sonwore   SONWORE      PPL
     adm1    lat_y    lon_x ave_den_allbuildings
1  BJ-000 11.24685 2.414291             804.0000
2  BJ-000 11.28007 2.427724             804.0000
3  BJ-000 11.28146 2.400951             804.0000
4  BJ-000 11.28761 2.381083             654.0000
5  BJ-000 11.28951 2.423543             804.0000
6  BJ-000 11.29832 2.401679             804.0000
7   BJ-AL 11.29845 2.438561             804.0000
8  BJ-000 11.30770 2.414829             595.8421
9   BJ-AL 11.33134 2.441141              92.1000
10 BJ-000 11.35827 2.422149              13.0000

10 Using points rather than non-intersecting buffers

Testing function for points:

path_points_gns <- file.path(path_to_project,
                        "INPUT","QGIS_map_voronoi","villages_points_adm4_gns.shp")

points_gns <- st_as_sf(st_read(path_points_gns))

points_gns <- st_transform(points_gns, 4326) #defining as EPSG4326 or WGS84
points_gns

Simple feature collection with 3599 features and 18 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 0.782811 ymin: 6.240975 xmax: 3.821388 ymax: 12.40508
Geodetic CRS:  WGS 84
First 10 features:
   ID gid_0   gid_1     gid_2       gid_3  name_1    name_2    name_3
1   1   BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
2   1   BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
3   1   BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
4   1   BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
5   1   BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
6   1   BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
7   1   BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
8   1   BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
9   1   BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
10  1   BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
           type_3 engtype_3 id_gns  full_name full_nm_nd sort_name desig_cd
1  Arrondissement   Borough   2058       Gano       Gano      GANO      PPL
2  Arrondissement   Borough   1426     Kokabo     Kokabo    KOKABO      PPL
3  Arrondissement   Borough    973      Wagui      Wagui     WAGUI      PPL
4  Arrondissement   Borough   3392   Atabénou   Atabenou  ATABENOU      PPL
5  Arrondissement   Borough   2726      Komon      Komon     KOMON      PPL
6  Arrondissement   Borough   2263      Dérou      Derou     DEROU      PPL
7  Arrondissement   Borough   3386  Banikoara  Banikoara BANIKOARA      PPL
8  Arrondissement   Borough    453 Toké-Banta Toke-Banta TOKEBANTA      PPL
9  Arrondissement   Borough   1873   Gomparou   Gomparou  GOMPAROU      PPL
10 Arrondissement   Borough    581    Sonworé    Sonwore   SONWORE      PPL
     adm1    lat_y    lon_x                  geometry
1  BJ-000 11.24685 2.414291 POINT (2.414291 11.24685)
2  BJ-000 11.28007 2.427724 POINT (2.427724 11.28007)
3  BJ-000 11.28146 2.400951 POINT (2.400951 11.28146)
4  BJ-000 11.28761 2.381083 POINT (2.381083 11.28761)
5  BJ-000 11.28951 2.423543 POINT (2.423543 11.28951)
6  BJ-000 11.29832 2.401679 POINT (2.401679 11.29832)
7   BJ-AL 11.29845 2.438561 POINT (2.438561 11.29845)
8  BJ-000 11.30770 2.414829  POINT (2.414829 11.3077)
9   BJ-AL 11.33134 2.441141 POINT (2.441141 11.33134)
10 BJ-000 11.35827 2.422149 POINT (2.422149 11.35827)

10.1 Calculating distance

ave_dist_rivers_points <- ave_dist_2shp(country_ext=gadm3,target=points_gns,shape_input=rivers_fao,UTM=31,res_raster=10000,name_dist="rivers")

[1] "UTM is: 31"

ave_dist_rivers_points[1:10, ] #print first 10 elements

   ID gid_0   gid_1     gid_2       gid_3  name_1    name_2    name_3
1   1   BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
2   1   BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
3   1   BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
4   1   BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
5   1   BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
6   1   BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
7   1   BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
8   1   BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
9   1   BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
10  1   BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
           type_3 engtype_3 id_gns  full_name full_nm_nd sort_name desig_cd
1  Arrondissement   Borough   2058       Gano       Gano      GANO      PPL
2  Arrondissement   Borough   1426     Kokabo     Kokabo    KOKABO      PPL
3  Arrondissement   Borough    973      Wagui      Wagui     WAGUI      PPL
4  Arrondissement   Borough   3392   Atabénou   Atabenou  ATABENOU      PPL
5  Arrondissement   Borough   2726      Komon      Komon     KOMON      PPL
6  Arrondissement   Borough   2263      Dérou      Derou     DEROU      PPL
7  Arrondissement   Borough   3386  Banikoara  Banikoara BANIKOARA      PPL
8  Arrondissement   Borough    453 Toké-Banta Toke-Banta TOKEBANTA      PPL
9  Arrondissement   Borough   1873   Gomparou   Gomparou  GOMPAROU      PPL
10 Arrondissement   Borough    581    Sonworé    Sonwore   SONWORE      PPL
     adm1    lat_y    lon_x ave_dist_rivers
1  BJ-000 11.24685 2.414291        5455.917
2  BJ-000 11.28007 2.427724        5483.610
3  BJ-000 11.28146 2.400951        4977.966
4  BJ-000 11.28761 2.381083        5037.752
5  BJ-000 11.28951 2.423543        5582.816
6  BJ-000 11.29832 2.401679        5329.104
7   BJ-AL 11.29845 2.438561        6046.705
8  BJ-000 11.30770 2.414829        5767.787
9   BJ-AL 11.33134 2.441141        6715.481
10 BJ-000 11.35827 2.422149        6887.261

10.1.1 Comparing results by polygons and by points

p1 <- hist(ave_dist_rivers_points$ave_dist_rivers)

p2 <- hist(ave_dist_rivers$ave_dist_rivers)

plot(p1, col=rgb(0,0,1,1/4))  # first histogram
plot(p2, col=rgb(1,0,0,1/4), add=T)  # add
legend("topright", legend=c( "ave_dist_rivers_points","ave_dist_rivers"), 
       fill=c(rgb(0,0,1,1/4), rgb(1,0,0,1/4)), border=NA)